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REACTIVE  SHOCK  PHENOMENA  IN  CONDENSED  MATERIALS: 
FORMULATION  OF  THE  PROBLEM  AND  METHOD  OF  SOLUTION 


I.  INTRODUCTION 

Shock  waves  are  particularly  useful  dynamic  phenomena  for  generating 
high  temperatures  and  pressures  very  quickly.  They  are  used  in  practical 
commercial  applications  such  as  synthesizing  diamonds,  in  military  applica¬ 
tions  such  as  igniting  explosives,  and  in  basic  research  applications  such  as 
producing  a  controlled  environment  for  studying  elementary  chemical  reactions 
rates.  In  this  paper,  we  describe  a  numerical  reactive  shock  model  developed 
at  NRL  which  is  being  used  to  study  the  physics  controlling  the  formation  and 
propagation  of  shocks  in  condensed  phase  materials.  When  coupled  to  a  chemi¬ 
cal  kinetics  reaction  scheme,  the  model  describes  the  detailed  chemical 
processes  initiated  by  shock  waves.  When  coupled  to  a  model  for  energy 
release,  the  code  describes  detonations. 

A  review  of  the  theory  and  equations  used  in  numerical  modeling  of 
shocks  and  detonations  in  condensed  phase  explosives  has  been  given  by  Mader 
[ 1 ] .  Mader  describes  the  physics  of  detonations  in  common  explosives  and 
presents  a  detailed  description  of  the  equations  of  state  of  the  condensed 
materials  and  gaseous  products,  the  model  we  have  developed  is  different  in 
several  ways  from  the  ones  he  has  developed  and  is  based  on  the  Reactive 
Shock  Models  described  by  Oran  and  Boris  [2].  The  covective  transport  equa¬ 
tions  are  solved  using  a  highly  accurate  finite-difference  algorithm,  FCT, 
which  has  both  minimal  diffusion  and  minimal  error  in  the  phase  propagation 
{3].  We  have  added  routines  for  external  energy  deposition  which  allow  us  to 
simulate  particular  mechanisms  such  as  laser  energy  input,  ihe  model  of  the 
fluid  dynamics  can  be  coupled  to  simple  reaction  models  or  to  detailed  chemi¬ 
cal  kinetics  subprograms  in  order  to  study  the  initiation  mechanisms  of  high 

energy  explosives.  Although  we  describe  a  one-dimensional  model,  the  basic 
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algorithms  for  the  equations  of  state,  energy  deposition,  and  chemical 
kinetics  can  be  easily  incorporated  in  two-  or  three-dimensional  FCT  codes 
and  so  the  calculations  are  directly  extendable  to  multi -dimensions. 

Finally,  the  methods  allcw  the  code  to  be  vectorized,  makinq  it  possible  to 
run  large  problems  efficiently  on  vector  computers. 

In  Section  II  we  describe  the  physics  and  chemistry  in  the  model  and  dis¬ 
cuss  the  formulation  of  the  problem.  First,  we  describe  the  equation  of 
state  used  in  the  model  for  condensed  phases.  An  Arrhenius  one-step  reaction 
supplemented  with  an  induction  model  is  assumed  to  control  the  conversion  of 
explosive  to  gas  products.  Then  we  describe  the  equation  of  state  for  the 
products.  When  both  phases  coexist,  temperature  and  pressure  equilibrium  are 
assumed.  In  Section  III  we  describe  the  numerical  solution.  Finally, 
typical  results  are  presented  in  Section  IV,  where  we  first  describe  a 
calculation  in  which  a  laser  pulse  deposits  energy  into  a  thin  layer  of 
absorbing  material  confined  between  a  wall  and  a  liquid.  A  shock  is  gener¬ 
ated  in  the  liquid  and  we  observe  its  characteristics  as  it  propagates  and 
decays.  We  then  describe  a  second  example  in  which  a  detonation  wave 
develops  from  a  hot  spot  in  liquid  nitrorae thane. 
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II.  FORMULATION  OF  THE  PROBLEM 


In  principle,  a  complete,  detailed  numerical  description  of  the  genera¬ 
tion  and  evolution  of  a  reactive  shock  in  a  homogeneous  material  involves  the 
solution  of  the  time-dependent  Navier-Stokes  equations  for  the  reactive  mate¬ 
rial.  We  also  must  incorporate  a  great  deal  of  information  about  the 
chemical  and  physical  properties  of  the  particular  material  involved.  We 
need  to  specify  the  equations  of  state  for  reactants  and  products  that  go 
properly  to  the  cold  solid  and  hot  gas  limits.  Heat  transfer  by  conduction 
is  usually  too  slow  to  affect  the  propagation  but  radiation  from  the  high 
temperature  products  to  the  undecomposed  reactive  material  surface  can  be 
quite  significant,  especially  when  the  products  contain  solid  carbon 
particles .  Thus  we  need  to  know  the  radiative  properties  of  the  gases 
produced  and  of  the  material  surface.  Also,  we  need  to  know  the  chemistry  of 
reactions,  or  at  least  the  rate  of  combustion  and  the  energy  release. 

Finally,  we  might  need  to  describe  other  effects  influencing  the  propagation 
process  such  as  species  diffusion,  carbon  deposition  and  coagulation. 

If  the  material  through  which  the  shock  propagates  is  heterogeneous, 
there  are  additional  complexities.  Most  commercial  explosives,  for  example, 
are  formed  by  compressing  the  explosive  powder  with  a  soft  binding  agent. 

When  such  an  explosive  is  subjected  to  stress,  its  deformation  properties  are 
complicated  by  the  yield  of  the  binding  agent.  Heat  conduction  becomes  an 
important  mechanism  for  the  growth  of  the  hot  spots  formed  by  wave  reflec¬ 
tions  at  the  complex  interfaces  between  the  powder  and  the  binding  agent.  It 
would  be  an  enormous  task  to  develop  and  use  a  model  that  would  treat  each 
element  of  an  explosive  powder  or  bindinq  material  as  a  separate  homogeneous 
domain  with  the  proper  compatibility  conditions  between  one  element  and 
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another.  Usually  this  problem  is  dealt  with  by  using  a  phenomenological 
composite  equation  of  state  whose  parameters  are  derived  experimentally  from 
measurements  of  shock  or  detonation  waves  propagating  in  the  actual  reactive 
material . 

In  this  work  we  consider  one-dimensional  shocks  in  homogeneous 
materials.  The  effects  of  bulk  viscosity  are  neglected,  therefore  omitting 
the  need  for  deviatoric  stress  terms  (since  shear  deformations  are  absent  in 
one  dimension).  Radiation  and  heat  conduction  are  assumed  insignificant  and 
therefore  neglected.  Therefore,  the  propagation  of  the  reactive  shock  is 
reduced  to  solving  the  set  of  equations  for  mass,  momentum  and  energy 
conservation  combined  with  the  proper  equation  of  state  and  a  description  of 
the  chemistry! 

( 1 )  Mass  Conservation 

3  .  <x—  ? .  3  ,  o—l .  .  . 

( pr  )  +  -gp  ( pu r  )  »  0  (i) 

( 2)  Momentum  Conservation 

3  ,  3—  1  3  7  3—  1  o-l  3p 

( pur  )  +  -gj-  (  pu  *r  )  +  r  -g£  «  0  (2) 

(3)  Energy  Conservation 

3  .  o—1  3  ,  o-l  3  o-l  o—1  • 

(Er  )  +  -g^  (Eur  )  +  -gj-  (r  pu)  »  r  pAl  (3) 

where  t  denotes  the  time,  r  the  radius,  p  the  density,  u  the  particle  veloc¬ 
ity  and  p  the  pressure.  The  total  energy  per  unit  volume  is  E  =  pl+pu2/2, 
where  I  is  the  specific  internal  energy  (including  the  stored  chemical 


energy)  and  Af  denotes  the  energy  absorbed  per  unit  mass  per  unit 
time  by  the  bulk  of  fluid  from  an  external  source,  the  exponent  a  is  a  geo¬ 
metric  index:  a  >  1,2  and  3  for  planar,  cylindrical  and  shperical  geometries, 
respectively. 

the  above  relations  are  supplemented  by  both  a  mechanical  equation  of 
state  and  a  caloric  equation  of  state,  written  symbolically  as 

P  -  P( P.I/W) 

T  ■  T( p,I,w), 

where  T  is  the  temperature  and  w  is  the  explosive  mass  fraction.  We  use  an 
equilibrium  equation  of  state  for  the  condensed  explosive  and  for  the  gaseous 
products .  thus  we  assume  that  the  different  intermolecular  and  intra¬ 
molecular  degrees  of  freedom  of  a  molecule  are  in  equilibrium  at  every 
instant  in  time.  In  other  words,  if  the  translational  energy  increases,  the 
rotational,  vibrational  and  electronic  degrees  of  freedom  absorb  enough  of 
the  new  energy  to  reach  the  same  temperature  instantaneously.  Althouqh  this 
assumption  is  not  strictly  correct  [4],  it  is  a  good  approximation  for  most 
purposes  and  it  enables  us  to  avoid  the  complexities  involved  in  a  model 
which  takes  the  relaxation  of  the  molecular  degrees  of  freedom  into  consider¬ 
ation. 

The  equation  of  state  for  the  condensed  explosive  is  usually  derived 
experimentally  by  determining  the  states  behind  explosion-generated  shocks  in 
the  explosive  material.  Since  the  pressure  generated  behind  the  detonation 
wave  in  the  condensed  phase  is  of  the  order  of  100  Kbars,  the  ideal  gas 
assumption  is  not  valid  for  the  gaseous  products.  As  a  conaquence,  we  cannot 
assume  that  the  different  gases  constituting  the  products  can  be  treated  as 
independent,  noninteracting  species.  A  theoretial  description  of  the 
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products  then  has  to  consider  the  forces  between  molecules  of  different 
species  as  well  as  those  between  identical  ones,  thus  adding  more  complica¬ 
tions  to  finding  an  equation  of  state  for  the  mixture  of  products.  The 
BXW  [5],  LJD  (6],  and  MCA  [7]  equations  of  state  are  all  attempts  to  describe 
mixtures  of  gasses  at  very  high  pressures . 

Finally,  we  need  to  know  the  rate  of  decomposition  of  the  material  in 
terms  of  the  instantaneous  state  parameters  such  as  temperature,  and  either 
pressure  or  density.  Specifically,  we  need  an  expression  for  the  evolution 
of  the  composition  with  time. 

II. 1.  Equations  of  State  of  the  Condensed  Phase 

The  first  equations  of  state  developed  for  condensed  phases  were  derived 
by  determining  the  isothermal  compressibility  using  a  method  developed  by 
Bridgman  [81 .  The  states  off  the  isotherm  are  then  calculated  from  thermo¬ 
dynamic  relations.  However,  practical  considerations  limit  isothermal  or 
static  compression  experiments  to  pressures  which  are  much  lower  than  those 
encountered  in  3hock  compressed  condensed  phase  materials,  especially  those 
occuring  in  detonations . 

In  response  to  this  problem,  Walsh  and  Christian  [9]  developed  a  tech¬ 
nique  for  deriving  an  equation  of  state  by  determining  the  locus  of  single- 
shocked  states,  i.e.  the  Hugoniot  curve.  States  off  the  Hugoniot  can  be  then 
evaluated  using  thermodynamic  relations.  The  Walsh  and  Christian  method 
involves  measuring  the  relation  between  the  shock  velocity  (produced  by  an 
explosion)  and  the  free  surface  velocity  in  the  explosive  material,  and  then 
using  the  fact  that  the  free  surface  velocity  is  almost  exactly  double  the 
particle  velocity  behind  the  shock  wave.  This  assumption  is  valid  to  1%  for 
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pressures  up  to  500  Kbar.  The  relation  between  the  shock  velocity,  ug,  and 
the  particle  velocity.  Up,  has  been  found  to  be  very  nearly  a  linear  fit 
[10] 


u  ■  c  +  s  u 
s  o  op 


(4) 


where  c_  is  the  sound  speed  in  the  undisturbed  material.  Once  c  and 
o  o 

sQ  are  measured,  the  pressure  of  the  condensed  phase,  denoted  from  now  on 
by  a  subscript  c,  is  given  by 


where 


*7  {I=  -  V 

c 


p  c2n 
o  o 

p  + - 

(i-s  n)2 


is  the  Hugoniot  pressure  corresponding  to  a  compression. 


(5) 


(6) 


n  = 


V  -  V 
o  c 

V 

o 


=»  ^ 


In  the  above  relations,  the  subscript  "o"  denotes  the  state  of  the  undis¬ 
turbed  material,  while  v  is  the  specific  volume  and  Y  is  the  Gruneisen  gamma. 
Hie  Hugoniot  internal  energy,  IH,  is  obtained  from 


Xo  +  2 


(P«  +  PQ)von. 


(7) 


Hie  temperature  Tc  is  given  by 
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T 

c 


(8) 


i  -r 

th 


where  TH  is  obtained  from 


In  T  =  P  +  G  (In  v  )  +  H  (In  v  ) 2  +  I  (In  v  )  3  +  J  (In  v  ) 4.  (9) 

Hsscsc  sc  sc 


The  coefficients  Fg,  Gg/  Hg,  I  and  Jg  of  the  expansion  are  derived  from  a 

least  squares  fit  to  the  Hugoniot  temperature  in  terras  of  the  logarithm 

of  the  volume.  The  details  of  the  derivation  of  the  condensed  phase  equation 

of  state  are  given  in  the  Appendix  A1 . 

Equations  (4)  to  (9)  are  used  when  v  <  v  ,  i.e.  for  compressed  states. 

c  o 

When  v  >  v  , 

c  o 


c 


(10) 


and 


T  *  T  +  ^^O  (11) 

c  o  — - 

V 

where  Cv  is  the  specific  heat  at  constant  volume  and  a  is  the  coefficient 
of  linear  expansion. 


II. 2.  Equation  of  State  of  the  Gas  Products  and  the  Rate  of  Reaction 

The  equation  of  state  of  the  gas  products  expresses  the  pressure  and 
temperature  in  terras  of  the  specific  volume,  the  internal  energy,  and  the 
composition  of  the  products.  If  the  detailed  chemistry  of  the  reaction  zone 
i3  sought,  the  reaction  zone  composition  could,  in  principal,  be  obtained 
theoretically  by  first  determining  and  then  solving  an  appropriate  set  of 
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chemical  kinetic  and  thermodynamic  equations.  On  the  other  hand,  if  we 
assume  an  infinitely  fast  rate  of  reaction,  the  reaction  zone  reduces  to  a 
discontinuity  propagating  through  the  explosive.  The  composition  of  the 
products  immediately  behind  the  wave  front  can  then  be  obtained  by  simulta¬ 
neously  solving  the  chemical  equilibrium  relations  and  the  jump  conditions. 
The  gas  products  then  expand  or  contract  isentropically  (maintaining  chemical 
equilibrium  at  every  point)  according  to  the  upstream  boundary  conditions. 

Since  very  little  is  known  about  the  mechanism  of  decomposition  of  con¬ 
densed  explosives,  our  preliminary  model  adopts  the  following  assumptions: 

a)  The  condensed  explosive  is  converted  to  qas  products  in  chemical 
equilibrium  according  to  the  global  finite  racticn  rate: 

dw  -E*/RT 

—  *  -  wZ*e  (12) 

b)  To  avoid  solving  the  chemical  equilibrium  relations  for  the  composi¬ 
tion  every  time  the  equation  of  state  of  the  products  is  needed,  we  assume 
that  the  states  of  the  gas  products  are  close  to  those  encountered  in  a 
steady  Taylor-type  detonation  wave.  In  this  type  of  wave,  the  condensed 
explosive  reaches  the  Chapman -Jouguet  (CJ)  equilibrium  state  instantaneously 
after  crossing  the  wave  front  and  then  expands  isentropically  to  the  undis¬ 
turbed  pressure.  The  equilibrium  isentrope  is  therefore  a  good  reference, 
the  same  way  the  Hugoniot  was  a  reference  for  the  compressed  condensed  explo¬ 
sive.  Then  we  can  use  an  averace  Gruneisen  gamma  to  evaluate  the  pressure 
and  temperature  of  states  off  the  isentrope.  The  details  of  the  derivation 
of  the  gas  products  equation  of  3tate  are  given  in  Appendix  A2. 
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The  pressure  of  the  gas  products,  denoted  by  a  subscript  g,  is  given  by 


p  *  p.  +  —  (I  -  I.  ) 
g  i  v  g  i 

g 


(13) 


The  pressure,  p^,  and  internal  energy,  1^ ,  on  the  isentrope  are  deter¬ 
mined  from 


In  p.  =  A  +  B( in  v  )  +  C(ln  v  )2  +  D(lnv  )  3  +  E(ln  v  ) 4  (14) 

l  g  g  q  q 


and 


ln(Ii+Z)  =  K  +  L(ln  pi)  +  M(ln  pj2  +  N(ln  p i)3  +  0(ln  p.^4.  (15) 


The  coefficients  of  the  expansions  in  Eg.  (14)  and  Eg.  (15)  are  derived  from 


a  least  sguares  fit  to  the  egui librium  isentrope.  in  the  above,  y?  is  the 


Gruneisen  gamma  evaluated  at  the  state  "i"  on  the  isentrope,  and  Z  is  a  con¬ 
stant  which  adjusts  the  standard  state  energy  to  be  consistent  with  that  of 
the  condensed  explosive.  Finally, 


T  =  T. 


I  -  I. 
+  1 


g  i  C* 
v 


(16) 


where 


In  T.  »  Q  +  R(ln  v  )  +  S (lnv  ) 2  +  T(ln  v  ) 3  +  U(ln  v  ) 4,  (17) 

1  g  g  g  g 


where  C'  is  the  average  specific  heat  of  the  products. 


II. 3.  Eguation  of  State  for  Mixtures  of  Condensed  Explosive  and  Gas 
Products 

Since  the  internal  energies  of  both  the  condensed  explosive  and  gas 
products  are  measured  from  the  same  reference  energy,  there  is  no  need  for  a 


heat  release  term  explicity  in  the  energy  equation.  Instead,  the  heat 


released  by  the  chemical  decomposition  appears  as  an  increase  in  temperature 
and  pressure.  In  other  words,  for  the  same  values  of  v  and  I,  the  equation 
of  3tate  of  the  products  yields  a  much  higher  temperature  and  pressure  than 
the  equation  of  state  for  the  condensed  explosive .  In  between  the  two 
extremes,  w  *  1,  representing  the  pure  condensed  explosive,  and  w  =  0,  repre¬ 
senting  the  gaseous  products,  the  temperature  and  pressure  increases  contin¬ 
uously  as  w  decreases. 

The  finite  rate  of  burning  expressed  in  Eq.  (12)  yields  a  reaction  zone 
of  finite  thickness  in  which  0  <  w  <  i  and  which  contains  a  mixture  of  con¬ 
densed  explosive  and  gas  products .  We  assume  that  in  those  regions  con¬ 
taining  such  a  mixture  the  condensed  phase  is  composed  of  a  large  collection 
of  small  fragments  uniformly  distributed  over  the  region.  Moreover,  we 
neglect  all  transient  effects  between  the  condensed  and  gas  phases  and 
assume  temperature  and  pressure  equilibrium,  Tc  *  and  Pc  *  Pg. 

Since  the  volumes  and  internal  energies  of  the  two  phases  are  additive, 
determining  the  equation  of  state  of  the  mixture, 

T  *  T(v,I,w) 
p  =»  p(v,I,w), 

reduces  to  finding  the  values  of  Ic,  1^,  vc  and  v^  that  satisfy 


T  (V  ,1  ) 

»  T  (v  ,1  ) 

(18a) 

C  C  C 

g  g  g 

p  (V  ,1  ) 

«  p  (v  ,1  ) 

(18b) 

c  c  c 

g  g  g 

I  *  wl  + 

(1-w)I 

(18c) 

c 

g 

V  *  wv  +• 

c 

(1-w)v 

g 

(18d) 

Equations  (18)  are  four  equations  in  four  unknowns  vc,  v^,  Ic  and 
1^.  They  can  be  reduced  to  a  single  nonlinear  equation  expressing  the 
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difference  in  pressure  (pc  -  p^)  in  terms  of  vc  (or  v^) .  To  obtain 

this  reduction,  we  first  enforce  the  temperature  equilibrium 

c_  *  T_  ■  T. 
c  g 


Then  we  rewrite  Eqs.  (8)  and  (16)  as 


Ic  =*  Ir(vc)  +  CV(T  -  Tr(vc)  ) 

i 

I  *  I. (v  )  +  C  (T  -  T  (v  )) 

Xg  2.  g  v'  i  g 

and  multiply  the  first  by  w  and  the  second  by  (1-w).  After  using  Eq.  (18c) 
we  solve  for  T: 


T 


I  +  w(CvTr  -  Ir)  +  (1  -  w)(C’T.  -  Ii) 

—  _w) 

V  V 


(19) 


In  Eq.  (19),  the  subscript  r  denotes  the  reference  state  of  the  condensed 

material.  As  explained  in  Appendix  A1 ,  r  denotes  H  when  n  >  0,  and  c' 

(v  3  v  ,  p  ,  =0)  when  n  <  0. 
c  c  c 

Next,  we  solve  for  the  pressure  difference  (pc-pg) •  Now  that  the 

equilibrium  temperature  is  known  for  a  given  choice  of  vc  (or  v^), 

substituting  for  (I  -I  )  by  C  (T-T  )  in  Eq.  (21)  and  for  (I  -I.)  by  C'(T-T. ) 

c  r  v  r  9  *  vi 

in  Eq.  (27),  we  get 


P  -  P„ 
c  g 


[Pr(vc)  -  Pi(vg)l  + 


YC  T?C- 

—  [T-T  (v  )]  -  -=—  [T-T.  ( v  )  ] 
v  r  c  v 

c  g 


(20) 


Equation  (2)  is  solved  by  iterating  on  vQ  (or  v^) . 


12 


II. 4.  Induction  Model 


Since  an  Arrhenius  one-step  reaction  alone  cannot  represent  an  extended 
period  with  no  apparant  temperature  or  pressure  change,  the  reaction  rate 
given  in  0g.  (12)  must  be  supplemented  with  a  model  for  such  an  induction 
time.  Thus  if  it  is  required,  an  induction  time  can  elapse  before  the  reac¬ 
tion  rate  of  Bq .  (12)  an  thus  the  energy  release  and  conversion  to  products 
is  switched  on.  If  the  steady  state  induction  time  (at  constant  temperature, 
T°,  and  pressure,  p°)  is  denoted  by  T°(T°,p°),  the  quasi-steady 
induction  time  for  varying  temperature  and  pressure,  t,  is  determined  from 
the  solution  of  the  integral  equation 


In  our  algorithm,  T  is  obtained  from  the  solution  of  f ( r)  =  1,  where 


df 

dt 


1 

T°(T,p) 


(21) 


and  f ( 0)  *  0.  Moreover,  T°  may  be  determined  emperically  or  using  detailed 
chemical  kinetics  calculations.  Here  we  will  assume  it  takes  a  simple 
Arrenius  forms 

,  £  /  o 

O  *  E  'RT  . 
x  ■  A  e 
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III.  NUMERICAL  SOLUTION 


We  must  solve  the  set  of  conservation  equations  (1)  to  (3),  the  rate  of 
reaction,  Eq.  (12),  and  the  induction  time  lapse,  Eq.  (21).  These  are 
supplemented  by  the  equations  of  state  for  the  range  0  <  w  <  1 .  First  we 
write  Eqs.  (12)  and  (21)  in  the  same  transport  form  as  Eqs.  (1)  to  (3).  This 
is  accomplished  by  transporting  pw  and  pf  instead  of  w  and  f,  namely 


3  ,  o-1 , 

lt(pwr  5 


3  o-l 

+  -^(  pwur  ) 


o-1  -E*/RT 

pr  wZ*  e 


(22) 


and 


3  .  ,  a-1 
*St<pfr  ) 


+  Pfur 


pr 


a-i 


T°  (p,  T) 


(23) 


Equations  (12)  and  (22)  are  equivalent,  as  are  Fqs.  (21)  and  (23).  The 
integration  step  starts  by  evaluating  the  time  step  5t  using  the  Courant 
condi tion 


6t  <  min  h — i — ) 

(uT+c 

where  c  is  the  sound  speed.  This  guarantees  the  stability  of  the  numerical 
solution  of  Eqs.  (1)-(3)  and  (22)-(23).  The  expressions  for  the  sound  speed 
in  the  condensed  phase,  gas  phase,  and  the  mixture  of  both  phases  are  given 
in  Appendix  A3. 
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the  solution  of  Eqs.  ( 1 ) — ( 3 )  and  (22)-(23)  is  obtained  by  operator 
splitting  the  fluid  dynamic  step,  the  chemical  reaction  step,  and  the  energy 
deposition  step.  'Bie  fluid  dynamic  step  is  solved  using  JPBFCT  (a  modified 
version  of  ETBFCT  [3]),  a  fourth  order  flux -corrected  transport  algorithm. 
Equations  (1),  (2),  and  the  fluid  dynamic  components  of  Eqs.  (3),  (22)  and 
(23),  the  left  hand  sides,  are  transported  using  JPBFCT  to  advance  p,  u,  I, 
w,  and  f.  The  explosive  mass  fraction  w  is  then  limited  between  0  and  1 
using 


w  =  max  {0,  min(1,w)} 


(24) 


and  the  equation  of  state  is  used  to  get  p  and  T. 

Next,  the  program  checks  if  the  induction  time  has  elapsed  (f  >1)  and 
if  so,  switches  on  the  chemical  reaction.  The  fraction  of  induction  time,  f, 
is  advanced  and  if  switched  on,  the  mass  fraction,  w,  is  also  advanced,  by 
solving  the  ordinary  differential  equations: 


3  ,  -  ohI 

g-  (  pfr  ) 


pr 


o-1 


T°(T,p) 


(25a) 


and 


3  .  o^1 

■St  (pwr  > 


0-1  * 
pwr  Z  e 


■E*/RT 


(25b) 


the  right  hand  sides  of  Eqs.  (22)  and  (23).  The  solution  for  f  is  given  by 
the  explicit  formula 
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I 


at 

T°(T,p) 


(26a) 


^new 


fold 


and  the  new  value  of  w  comes  from  the  implicit  formula 


new  old.,.  *  -E  /RTV, 
w  *  w  /  ( 1 +Z  e  5t ) , 


(26b) 


Hie  equation  of  state  is  then  used  again  to  determine  p  and  T,  after  the 
burning  process. 


Finally,  if  energy  is  deposited  from  an  external  source,  it  is  added  to 


the  specific  internal  energy  I,  and  the  equation  of  state  is  used  again  to 
obtain  p  and  T  in  preparation  for  a  new  integration  time  step. 
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IV.  TYPICAL  RESULTS 


IV. 1.  Laser-Induced  Shock  Wave  Structure 

The  calculation  presented  below  was  performed  to  assess  the  ability  of 
long  duration  laser  pulses  to  create  a  shock  wave  structure  with  a  reasonably 
wide  plateau  behind  the  shock.  Since  the  plateau  would  provide  an 
environment  with  constant  pressure  and  temperature  in  which  optical  probes 
measurements  would  be  relatively  easy  to  analyze,  it  is  a  potentially  impor¬ 
tant  feature  of  an  experiment.  In  an  upcoming  report  we  will  discuss  those 

properties  of  the  material  and  laser  pulse  which  control  the  evolution  of  the 
plateau.  In  this  section  we  describe  one  test  calculation  in  some  detail. 

Consider  a  30ym  layer  of  plexiglass  confined  on  one  side  by  a  rigid  wall 
and  on  the  other  by  a  semi -infinite  slab  of  water.  A  20  ns  half -width  laser 
pulse  is  deposited  into  the  plexiglass  in  such  a  way  that  the  energy  depos¬ 
ited  per  unit  volume  decays  linearly  with  the  depth  into  the  plexiglass  and 
vanishes  at  the  edge  where  the  plexiglass  meets  the  water. 

The  computational  grid  spacing  in  each  material  is  uniform,  but  varies 
with  the  material.  The  initial  cell  size  in  the  plexiglass  zone  is  taken  as 
1 .5  urn.  The  initial  cell  size  in  the  water  zone  is  then  determined  such  that 

the  ratio  of  the  cell  sizes  is  nearly  equal  to  the  ratio  of  the  local  sound 

speeds.  Thus  the  Courant  time-step  is  nearly  equal  in  both  regions.  The 
number  of  cells  in  each  zone  is  fixed.  The  interface  of  the  last  cell  in  the 
plexiglass  zone  moves  with  the  local  particle  velocity.  As  a  result  the 
cells  in  the  plexiglass  zone  expands  while  those  in  the  water  contract. 
However,  each  region  maintains  uniform  grid  spacing  at  all  times.  The  power 
of  the  laser  pulse  is  assumed  to  be  triangular  in  shape  with  a  40  ns  base 
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width.  The  initial  temperature  and  pressure  are  300°K  and  1  atm, 
respectively.  The  parameters  used  in  the  different  equations  of  state  are 
given  in  *ftible  I. 

The  calculated  pressure  at  different  times  is  illustrated  in  Fig.  (1). 
The  position  of  the  interface  between  the  plexiglass  an  water  is  denoted  by  a 
dot  on  each  of  the  shown  profiles.  As  the  input  power  increases  for  the 
first  20  ns,  the  pressure  near  the  wall  increases  continuously.  When  the 
power  begins  to  decrease,  the  pressure  at  the  wall  starts  to  drop  and  a 
propagating  wave  develops  (26  ns).  The  calculation  predicts  a  10-20  ns 
plateau  behind  the  shock  wave  with  a  peak  pressure  which  is  essentially 
constant  for  a  time  longer  than  the  duration  of  the  laser  pulse  (44  ns) . 

Later  in  the  calculation,  the  details  of  the  initial  energy  deposition 
process  are  essentially  forgotten.  Then  the  wave  continously  approaches  the 
classical  blast  wave  profile  which  would  result  if  a  fixed  amount  of  energy 
is  deposited  instantaneously  in  an  infinitismally  thin  layer  of  water.  Such 
a  blast  wave  calculation  is  illustrated  in  Fig.  (2).  Here  we  notice  that  the 
shock  wave  pressure  decays  faster  than  in  Fig.  ( 1 ) . 

The  plateau  results  from  the  interaction  betwen  the  way  the  laser  energy 
is  deposited  and  the  way  it  is  transported  to  the  water.  In  Fig.  (1),  we 
notice  that  during  the  power  rise-time  of  the  laser  pulse  and  up  to  22  ns, 
the  pressure  profile  decays  away  from  the  wall.  During  this  period,  the 
pressure  builds  up  in  the  plexiglass  faster  than  the  material  can  expand. 
While  the  pressure  is  building  up  within  the  plexigless,  compression  pressure 
waves  emanate  from  this  region  and  propagate  into  the  water.  As  a  result, 
the  pressure  profile  follows  the  trend  of  the  energy  deposition  curve. 
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TABLE  I 


Plexialass 

Water 

p  *1.18  gm/cc 
o 

p  *  1 .0  gm/cc 
o 

c  *  0.243  cm/ps 
o 

c  *  .01483  cm/ us 

o 

s  -  1 .5785 
o 

s  -1.97 
o 

P 

M 

o 

1 

•r 

o 

n 

*— • 

o  »  6  *  10  J  K  1 

Y  *  2.157 

Y  *  1 .65 

o 

C  *  0.35  cal/gm  K 

C  *1.0  cal/gm  °K 
V 

v 
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p  (Kbar) 


r  (MICRONS) 

Figure  1  Calculated  pressure  profiles  for  a  20  ns  half -width  laser  pulse 

deposited  in  30  microns  of  plexiglass  in  contact  with  water.  Hie 
dot  on  each  curve  denotes  the  position  of  the  interface  between 
plexiglass  and  water. 
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p  (Kbar) 


r  (MICRONS) 

Figure  2  Classical  blast  wave  solution  in  water  for  the  same  amount  of 
energy  deposition  as  in  Fig.  1. 
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When  the  power  begins  to  drop,  the  pressure  in  the  plexiglass  drops  even 


faster  due  to  material  expansion.  A  rarefaction  wave  is  created  near  the 
wall  and  propagates  into  the  water,  while  the  compression  waves  downstream 
steepen  up  to  form  a  shock.  The  rarefaction  and  the  compression  waves 
combine  to  create  a  zone  of  relatively  uniform  pressure  right  behind  the 
shock.  During  this  period  the  peak  pressure  remains  nearly  constant  for  a 
period  of  time.  Eventually  at  40  ns  in  Fig.  (1),  the  rarefaction  wave  catch 
up  with  the  shock  resulting  in  a  plateau  immediately  behind  it,  and  later 
(51  ns)  the  plateau  disappears.  As  the  shock  wave  propagates  downstream,  the 
rarefaction  waves  cause  the  profile  to  change  from  convex  to  concave,  and 
approach  the  classical  blast  wave  profile. 

Finally,  Fig.  (3)  illustrates  the  calculated  pressure  waves  for  the  same 
conditions  of  Fiq.  (1),  except  that  all  cell  sizes  are  twice  as  large.  This 
serves  as  a  test  of  the  convergence  of  the  numerical  solution.  A  close 
comparison  between  the  two  figures  reveals  that  the  pressure  waves  corre¬ 
sponding  to  the  same  time  are  nearly  coincident.  Those  of  Fin.  (3)  are 
slightly  lower  in  value  near  large  pressure  gradients  due  to  the  larger  nume¬ 
rical  diffusion  associated  with  laroer  cell  sizes. 

IV. 2.  Evolution  of  a  Detonation  Wave  from  a  Hot  Spot 

The  next  calculation  shows  that  heating  a  high  energy  explosive 
uniformly  does  not  necessarily  lead  to  homogeneous  ignition.  Ignition  may 
start  at  hot  spots  caused  by  the  inhomogenieties  of  either  the  material  or 
the  heating  mechanism.  The  reaction  would  then  propaaate,  consuming  the 
material  between  hot  spots  before  any  significant  reaction  occurs  in  the 
background  material  which  is  at  slightly  lower  temperature. 
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p  (Kbar) 


Figure  3 


0  50  100  ISO  200  250 

r  (MICRONS) 

Calculated  pressure  profiles  for  the  same  physical  conditions  as 
in  Fig.  1,  but  for  computational  cell  sizes  twice  as  large. 
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Here  we  model  a  sample  of  liquid  nitromethane  which  is  1.2  cm  long.  It 
is  assumed  that  the  sample  is  initially  at  1000°K.  This  is  beyond  its 
iqnition  temperature  so  that,  in  principle,  it  should  start  reacting.  Hie 
pressure  is  initially  1  atm  yielding  a  background  density  of  0.692  gm/cc. 

Hie  hot  spot  is  modelled  by  superimposing  a  Gaussian  temperature  profile 
200°  K  higher  than  ambient  and  0.1  cm  wide  on  the  background  temperature  at 
the  center  of  the  sample.  This  symmetry  allows  us  to  perform  the 
calculations  for  half  the  sample  only,  and  the  center  of  the  Gaussian  is 
formally  equivalent  to  a  rigid  wall. 

Hie  computational  grid  is  uniformly  spaced  with  a  fixed  cell  sice  of 
0.01  cm.  Hie  parameters  of  the  eqution  of  state  and  the  rate  of  reaction  for 
liquid  nitromethane  are  given  in  Table  II.  We  notice  that  A*  =  0.  In  this 
case,  it  is  assumed  that  the  parameters  of  the  one-step  Arrhenius  reaction 
are  adjusted  to  crudely  incorporate  any  induction  time. 

Hie  calculation  predicts  the  formation  of  a  detonation  wave  within 
0.4  us.  As  illustrated  in  Fig.  (4),  the  detonation  pressure  builds  up 
quickly  towards  a  steady  wave  pressure  of  65  Kbars  that  propagates  at 
0.54  cm/ us  (Mach  number  =  3.72).  At  1.28  us,  the  detonation  wave  initiated 
at  the  hot  spot  has  almost  consumed  the  whole  sample  before  any  significant 
reaction  occurs  at  the  background  temperature.  Without  the  presence  of  the 
hot  spot,  the  material  would  have  exploded  31.2  us,  much  later  than  it  has 
here. 
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TABLE  II 


Liquid  Nitromethane 


p  ==  1.1 28  qm/c 
o 

c  *  0*1  647  cm/  {is 
o 

s  =  1 .637 
o 

a  =  3  *  1 0’"1*  °K-1 
Y  =  0.6805 

=  0.414  cal/qm°K 
=  0.556  cal/qm°K 
3  =  0.1 

2*  =  4.0  x  108  us"1 
E*  =  5.36  x  104  cal/mol 


A  =  -3.11585 
B  =  -2.35968 
C  *  +2.10663  x  10"1 
D  =>  +3.80357  x  10"3 
E  -  -3.53455  x  10“3 


Fg  =•  +5.41171 
Gg  =  -2.72959 
Hg  =  -3.21986 
Ig  =  -3.90757 
Jg  a  +2.39028 

K  =  -1 .39937 
L  =  +4.79350  x  10"1 
M  »  +6.06708  x  10"2 
N  =  +4.10673  x  10"3 
0  a  +1 .13327  x  io"4 


q  =  +7.79645 
R  a  -5.33007  X  10"1 
S  -  +7.09020  x  10"2 
T  *  +2.06150  x  10"2 
U  »  -5.66140  x  io"3 
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APPENDICES 


A1 .  Equation  of  State  of  the  Condensed  Phase 

Here  we  adopt  the  Walsh  and  Christian  method  [9]  for  deriving  the 

equation  of  state  of  a  condensed  explosive.  Wiis  method  begins  by  fitting 

the  measured  shock  velocity,  u_,  and  measured  particle  velocity,  u  .  to 

a  P 

the  linear  relation, 

U  ■  C  +  3  U  . 

3  O  OP 

Then  combining  the  above  equation  with  statements  of  conservation  of  mass 


-  V°H 


u  p 

3  O 


and  conservation  of  momentum 


p_  -  P  *  p  u  u 
H  O  O  S  p 


across  the  shock,  we  get 


pc  2(  1  -  ~) 


o  o 


PH  -  Po  - 


H 


(1  -3  (1  -  -p)]2 

°  PH 


(A1  .1  ) 


where  the  subscript  "o"  denotes  the  state  of  the  undisturbed  material,  and 
the  subscript  "H"  denotes  a  point  on  the  Hugoniot.  Denoting  the  volumetric 
compression  by  n, 


n  = 


(A1.2) 
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where  v  is  the  specific  volume,  Bq .  (A1 .1 )  reduces  to 


PH 


p  c2n 
o  o 

—  e 

M-s  n)2 

O 


(A1 .3) 


Finally,  from  energy  conservation  across  the  shock  front,  the  internal 
energy,  IH,  can  be  written  as 


:H  -  Zo  "7  (PH  +  PC)(Vo  -  V  ml  'PH  +  P0>  V* 


(A1 .4) 


3y  using  the  fact  that  temperature  has  a  relatively  small  effect  on 
compressibility,  Walsh  and  Christian  [9]  were  able  to  determine  the  variation 
of  the  temperature  alone  the  HUgoniot  by  assuming  that  Cv  and  Op/3T)v 
are  constants! 


—  bv  rT 

_  _  b(v  -  V,J  H  r  H  bv 

T„  »  T  e  o  H  ♦  e  J  e  f(v)dv 

H  o  — — —  'v 

C  o 


(A1 .5) 


where  Cv  is  the  specific  heat  at  constant  volume.  In  Bo.  (A1.5), 


while 


ft  >  -  1  dP” 

"V  -73^ 


<Vo"'h>  *  7  PH 


b  = 


(3p/3T)v 

C 

v 


By  performing  the  integration  in  Bq.  (A1.5),  we  obtain  the  temperature  along 
the  HUgoniot.  Since  in  a  typical  calculation,  Bq.  (A1 .5)  has  to  be  evaluated 
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each  time  the  temperature  is  needed,  Mader  [1]  fit  the  TH  - 
using  a  fourth  degree  polynomial  in  log  of  the  volume: 


vR  curve 


in  T. 


G  (in 
s 


V 


H  (  in 
s 


V 


I  (  in 
3 


V 


J  ( in 
s 


V 


( A1  .6) 


aquations  (A1.1)  to  (A1.5)  are  quite  general  and  are  equally  valid  for 
3olids  and  liquids.  For  some  solids,  however,  a  kink  in  Rq.  (4)  occurs  at 
the  point  where  there  is  a  change  of  phase  to  liquid.  In  this  case,  another 
set  of  coefficients,  Cj  and  s^,  are  required  in  addition  to  cq  and  sq. 

In  reactive-shock  phenomena,  the  3tates  attained  after  compression  are 
expected  to  fall  near  the  Rankine-Hugoniot  locus,  making  this  line  a  good 
reference  state.  Thus  with  the  use  of  the  Gruneisen*  Y,  defined  as 
Y  =  v(3p/3l)v,  we  can  reach  any  state  off  the  Hugoniot  from  point  H  which 
has  the  same  specific  volume  as  this  state.  This  is  illustrated  schemati¬ 
cally  in  Fig.  (Al.a).  Denoting  the  condensed  explosive  by  the  subscript  "c", 
we  have  for  constant  Y, 


PC  ’  PH(V 


—  [I  -  I  (v  ) 
V  C  He 

c 


(A1 .7) 


Similarly,  since  C 


(21, 

3T  v' 


for  a  constant  C 

v 


T 

c 


-  vv 


I  - 

c 


W 


(AT .8) 


•For  an  ideal  gas  Gruneisen  gamma,  Y  5  v(2E)  »  ( ilSZ)  (s  denotes  the 

3l  v  3lnv  s 

entropy)  is  related  to  the  adiabatic  index  Y* a  =  (2i2E,  by  y  *  Yj„-1 . 

31nv  s 
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Figure  A1  (a)  Path  used  in  deriving  the  equation  of  state  of  a  compressed 
condensed  explosive  beginning  from  its  undisturbed  state. 

(b)  Path  used  in  deriving  the  equation  of  state  of  an  expanded 
condensed  explosive  beginning  from  its  undisturbed  state. 


where  pu(v  ),  Ia(v  ),  and  T„(v  )  are  evaluated  using  Bqs.  (A1.3),  (A1.4),  and 
(A1.6)  when  vH  ■  vc. 

Equations  (Al.i)  to  (A1 .8)  are  used  to  determine  the  state  of  a 
compressed  condensed  explosive.  For  an  expanded  explosive,  however,  a 
different  method  is  used  to  obtain  the  state  properties.  This  is  illustrated 
in  Fig.  (Al.b).  Since  going  from  o  to  o*  is  a  constant  volume  process,  we 
obtain 

1,-1  -  PoVo  ,  ( A1 .9) 

°  Y 

where  p'»  0  by  definition.  Also, 
o 

T  .  -  T  +  ro'  -  Xo  .  (A1  .10) 

°  °  — c - 

v 


Mew  the  path  o'  to  c'  is  a  constant  zero-pressure  process,  so  that  the 
enthalpy,  H,  is  equal  to  the  internal  energy,  1.  We  then  can  write 


OI/3T)  (3H/3T) 

t-JL)  m  P  m  .  ..  .  _ E 

' 3v  (3v/3t)  (3v/3t) 

P  P  P 


The  numerator  is  the  specific  heat  at  constant  pressure,  Cp,  while  the 

denominator  is  related  to  the  coefficient  of  volummetric  expansion 
1  3v 

a  =  — (-s=)  .  For  zero  pressure,  C  2  c  .  Moreover,  for  condensed  phases 
v  v  3T  p  e  pv 

(solid  or  liquid),  vc  is  expected  to  be  close  in  value  to  vQ.  Thus 


<£>_  * 


3v  p  "  3  a  v 


where  a  is  the  coefficient  of  linear  expansion  ( <*v  »  3a). 


(Ai .1 1 ) 
Integrating 
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( Al  .1 1 ) ,  we  obtain 


B*. 


I  .  +  V  (v  -  V  ) 

o'  — -  c  o 

3a  v 

o 


Also,  since  Cp  2  Cv  along  o' 


-  =' 


T  »  T  ,  + 
c*  o' 


I  ,  I  , 
e  •  -  o' 


(Al .1 2) 


(Al .1 3) 


Combining  Eas.  (A1.9)  and  (A1.12)  yields 


1 


c' 


P  v  c  n 

o  o  V  . 

o  "  Y  3a 


(Al .14) 


Combining  Bqs.  (A1.10)  and  ( Al .13)  yields 


T 

c 


I 


V 


Finally  for  the  constant  volume  process  c'  -  c, 


(Al  .15) 


and 


c 


(Al  .16) 


T  (v  )  + 
c'  c 


I  -  I  ,(v  ) 
c  c*  c 


(Al  .17) 


Equations  (A1.14)  to  (A1.17)  are  referred  to  as  the  Gruneisen  equation  of 
state . 
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In  summary,  p  (vc,Ic),  T(vc,ic)  are  expressed  in  terms  of  a 
reference  state,  denoted  from  now  on  by  a  subscript  r,  namely 


p( V  ,1  )  «  p  (v  )  +  —  (I  -  I  (V  ) 
c  a  rc  v  crc 
c 


(A1  .18) 


T(v  ,  I  )  *  T  (v  )  + 
c  c  r  c 


I  (v  ) 
r  c 


(A1 .19) 


For  v  <  vQ,  r  =  H;  while  for  v  >  vq  ,  r  =  c' . 

It  remains,  finally,  to  define  a  standard  state  at  which  the  internal 
energy  is  set  equal  to  zero  by  definition.  Here  we  pick  IQ  =  0.  As  men¬ 


tioned  earlier,  this  is  the  undisturbed  state,  namely  p  =  1  bar,  T  2  300K. 

o  o 


A2.  aquation  of  State  of  the  Gas  Products 

The  gas  products  are  assumed  to  reach  chemical  equilibrium  instantane¬ 
ously  after  their  production.  To  solve  for  an  equilibrium  state  we  need  an 
equation  of  state  for  the  mixture  of  aas  products  in  terms  of  the  composi¬ 
tion.  According  to  Cowan  and  Fickett  [11],  the  BKW  equation  of  state  for  a 
mixture  of  gases  can  be  written  in  the  form 


where 


V*  ,  3x 

RT  *  1  +  x  e 


(A2.1  ) 


k 

v(T+e)a 


k  being  the  average  covolume  defined  in  terms  of  the  individual  covolumes, 
k^,  as  k  =  x^k^,  where  xi  is  the  molar  fraction  of  species  i.  the  a,  0, 

9  and  <  are  constants  adjusted  to  reproduce  the  detonation  Chapman-Jouquet 
pressure  and  velocity  obtained  experimentally.  If  we  know  the  internal 
energy  and  entropy  of  the  mixture  of  products  at  the  ideal  aas  limit  (p=*0), 
we  can  extend  their  validity  to  the  hundreds  of  Kilobars  pressure  range  by 
using  the  thermodynamic  relations 


(— — )  =*  T(— -  p 

3v  T  T  3t'v  p 


,  3s.  ,3v\ 

3p  T  dT  p 


together  with  Bq.  (A2.1).  The  equilibrium  state  is  then  obtained  by 
minimizing  the  Gibbs  free  energy.  Note  that  because  of  the  possible  presence 
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of  solid  carbon  in  the  products,  a  solid  equation  of  state  for  graphite  is 
needed.  According  to  Cowan  and  Fickett,  the  carbon  pressure  is  expressed  in 


terms  of  a  second  degree  polynomial  in  T,  the  coefficients  of  which  are  them¬ 
selves  polynomials  in  the  compression  of  the  solid  carbon  relative  to  its 
normal  crystal  density. 

To  avoid  solving  the  chemical  equilibrium  relations  every  time  the  equa¬ 
tion  of  state  of  the  products  is  needed,  it  assumed  that  the  states  encoun- 
terd  in  a  typical  calculation  are  close  to  those  encountered  in  a  steady 
Taylor-type  detonation  wave,  i.e.  states  on  the  equilibrium  isentrope  through 
CJ  point.  Mader  [1]  fit  the  equilibrium  isentrope  throuah  the  CJ  point  for 
different  explosives  to  a  set  of  polynomials  of  the  form 

in  p.  =  A  +  B  (in  v  )  +  C(  in  v  )2  +  0(  £n  v  )3  +  E(  in  v  ) 4  (A2.2) 

l  q  a  q  q 

in  r  =  K  +  L(  in  pj  +  M(  in  p^2  +  M(  in  p^3  +  0(  in  p^4  (A2.3) 

in  T.  ■  5  +  R(in  v  )  +  S(!n  v  )2  +  Kin  v  )3  +  U(ta  v  )4.  (A2.4) 

1  g  g  g  g 

The  shifted  internal  energy  +  Z,  where  Z  is  a  constant  used  to 

change  the  standard  state  to  be  consistent  with  the  condensed  explosive  one. 

For  points  off  the  isentrope,  the  definition  of  Gruneisen  Y'  and  c^  (for 
gases  we  use  an  apostrophe)  provide  the  transition  from  point  i  to  g,  as 
shown  in  Fig.  (A2)  .  Since  we  need  the  equilibrium  ?(v,I)  and  T(v,I),  and  we 
do  not  need  to  obtain  the  composition  of  the  products,  the  change  in  composi¬ 
tion  along  the  constant  volume  process  i-g  in  Fig.  (A2)  can  be  absorbed  in 
the  definition  of  Y'  and  C^.  As  toz  V,  since  g  is  expected  to  be  near  i,  we 

can  take  Y'  *  Y.'  where  y*  takes  into  account  the  change  in  5p/5l  due 

average  i  i  ^ 
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to  change  in  composition  along  the  equilibrium  constant  volume  process. 
First,  we  write 


i  v(|£) 
31  v 


^v 


(4J> 

3s  v 


(A2.5) 


then  substitute  for  (3p/3s)v  and  (31/ 3s)  from  the  first  law  of 


thermodynamics  (dl  *  Tds  -  pdv),  yielding 


V 


i  3<I\ 

v(-  -r~ ) 
3v  s 


/T  =  ( 


3lnT. 
3lnv  s 


Thus  y!^  is  expressed  as 
dinT. 

Y!  =>  ■— f — -  *  -  {R  +  2S(  inv  )  +  3T(  fnv  )  2  +  4U(  JJnv  )  3}  (A2.6) 

l  dJlnv  g  g  g 

g 


therefore  inferring  the  change  in  compositon  alonq  the  constant  volume 
process  i-g  in  Fig.  (A4)  from  the  change  along  the  constant  entropy  process 
through  i,  expressed  in  Bgs.  (A2.4).  Integrating  Eq.  (A2.5),  we  get 


Yi 

W  +”  (Ig  -  VV 

g 


(A2.7) 


Finally,  we  assume  that  is  insensitive  to  the  change  in  composition, 
giving 


I  -  I.  (v  ) 

t  ■  t ,  ( v  )  +  -2 — — 2-  .  (A2.8) 

g  i  g  cr 
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A3 .  Expressions  for  the  Sound  Speed 


Generally,  the  sound  speed,  c,  is  defined  as 


-*2<£>s 


Since  p  p(v,I),  let  I  *  I(v,s),  so  that 


'll'!* 


(i£\  r-ii) 

Jl  »  3v  s 


Substituting  for  (3l/3v)s  by  -p  (dl  *  Tds  -pdv)  we  get 


c2  = 


v*<£) 

3v  s 


3p 


I?, 


l*(-£V  Wt1 


(A3 .1  ) 


Condensed  Phase 

ttie  pressure  is  expressed  in  terms  of  a  reference  pressure,  Pr,  and  a 
reference  internal  energy,  Ir,  both  functions  of  specific  volume  only: 

p  =  p  (v)  +  —  ( I  -  I  ( v) )  . 
r  v  r 

Substituting  the  above  relation  in  Eq.  (A3.1)  we  obtain 

dp  dl 

c2  =■  v[ YP  -  V  +  (p  -  Pr)  +  y  1  .  (A3. 2) 

Since  the  reference  state  r  depends  on  whether  the  condensed  explosive 
is  compressed  or  expanded  we  consider  the  two  cases: 
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(1)  v  <  vQ  ,  r  =  H:  Protn  Bq .  (A1 .3)  and  the  definition  of  n  in  Bq. 

( A1  .2) , 


dv  “  dv 


v 

o 


+  sn. 


-  sn 


(pc  ) 2( i  +s  n) 

o  o _  o 

(i-s  n)3 

O 


(A3. 3) 


whereas  Eas.  (A1.2)  and  (A1.4)  give 

dIr  dIH  1  dpH 

d^"2—  ~ffld7--  (PH+  Po)po]  (A3‘4> 

o 

At  the  undisturbed  state  (pQ,  TQ),  n  »  0.  Bq.  (A3. 2)  gives  then 
c2  *  cq2,  as  shouid  be.  Note  that  cq  and  sq  are  replaced  by  c[( 
when  dealing  with  a  solid  that  exhibits  a  kink  in  the  linear  relation  between 
the  particle  velocity  and  shock  velocity  at  the  point  where  a  phase  change  to 
a  liquid  occurs . 

dpr 

(2)  v  >  v  ,  r  =  c':  Since  p  ,  =  0,  - —  *  0.  From  Bq.  (A1.14), 
o  c  av 

dl  dl  ,  Cp 

r  c 1  v  o 

■3 —  =  j .  At  the  undisturbed  state  (p  ,  T  ) ,  Bq.  (A3. 2)  yields 
av  av  30i  o  o 

yc 

c2  «  (Y+1)  PoVq  +  ^  , 


thus  giving  an  estimate  of  cc  as 


c 

o 


(( Y+1 ) 


P  v 
o  o 


♦ 


_ 1/2 

3a  1 


(A3. 5) 


For  liquid  nitromethane,  for  example,  Mader  [i]  gives  the  following  set  of 
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parameters 


C 

v 


cQ  =  1.647  x  io-1  cm/usec 
Y  *  6.805  x  10"1 


3  x  IQ-4  °K~1 


.  ..  „  ..-1  cal  4.14  xio“l  ,  .? 

4.14  x  io  — - 23890  (c"/Usec) 

q  m  X 


YCv  1/2  , 

so  that  [ ( Y+1 )p  v  +  - -}  =  1.145  x  io~l  cm/usec.  The  discrepancy  between 

o  o  3a 

cq  given  and  that  evaluated  from  Eq.  (A3. 5)  indicates  a  fcink  in  the  thermo¬ 
dynamic  properties  at  the  interface  between  the  two  reqimes  v  <  vQ  and 
v  >  v0.  However,  since  the  sound  speed  is  only  used  to  derive  the  Courant 
time  step,  which  is  an  upper  limit  for  the  actual  time-step  used,  this 
discrepancy  presents  no  problem  in  the  calculations . 


Gas  Phase 

The  equilibrium  isentrope  through  the  CJ  point  is  the  locus  of  reference 
states.  Upon  substituting  Eq.  (A2.7): 

Y’(v) 

p  =  pi(v)  +  —  [I  -  Ii(v)] 

in  Bq,  ( A3 • 1 ) ,  we  get  a  relation  similar  to  Eq.  (A3. 2)  except  for  a  correc¬ 
tion  term  attributed  to  the  dependence  of  on  v,  namely 

dp.  dl.  dy* 

c2  „  v[y.p  -  +  (p  -  p.)  +  Y.  ]  .  ( A3 . 6 ) 
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Noting  that  Zn  and  Y?  are  given  in  terms  of  Jin  v,  whereas  Jin  I£  is  given 
in  terns  of  in  p^,  Eq.  {A3. 6)  can  be  rewritten  as 


vtY^P 


dinp^ 

dinv 


+  (P 


Pi,] 


Y'(I.  + 


2) 


dJlnl^  dJlnp^ 


dlnp.  dlnv 


(A3. 7) 


dinp.  dlnl!  dY.' 

where  — — ,  — -  ,  and  — : —  are  obtained  by  direct  differentiation  of 

dJlnv  dinpi  dinv 

Eqs .  (A2.2),  (A2.3),  and  (A2.6),  respectively. 


Mixture  of  Condensed  Phase  and  Gas  Phase 

Since  the  sound  speed  in  the  mixture  is  known  a  priori  to  be  between  the 
values  of  sound  speeds  of  the  condensed  phase  and  gas  products,  it  is 
satisfactory  to  assume 


c  .  _  *  max(c  ,c  ] 

mixture  c  g 


{A3 .9) 


when  determining  the  Courant  time  step.  Here,  c_  and  c  are  evaluated  at 
the  common  temperature  and  pressure  of  the  mixture. 
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